fish_traits_sum <- fish_traits %>%mutate(across(where(is.numeric), round, 2))## Display the table ----htmltools::tagList(DT::datatable(fish_traits_sum))
traits correlation
Code
M <-cor(numeric_traits[, c(-1)])ggcorrplot::ggcorrplot(M, hc.order =TRUE, type ="lower",lab =TRUE, tl.cex =9, lab_size =3)
Code
# list of species sp_names <-c(rownames(fish_traits))# taxonomic_familiestaxonomic_families <- sp_names %>%as.data.frame() %>%`colnames<-`("species") %>%mutate(family =case_when( species %in%c("Benthosema_glaciale","Ceratoscopelus_maderensis","Diaphus_metopoclampus","Lampanyctus_ater","Lampanyctus_crocodilus","Lampanyctus_macdonaldi","Lobianchia_gemellarii","Myctophum_punctatum","Notoscopelus_bolini","Notoscopelus_kroyeri","Bolinichthys_supralateralis" ) ~"Myctophidae", species %in%c("Borostomias_antarcticus","Chauliodus_sloani","Malacosteus_niger","Melanostomias_bartonbeani","Stomias_boa" ) ~"Stomiidae", species %in%c("Holtbyrnia_anomala","Holtbyrnia_macrops","Maulisia_argipalla","Maulisia_mauli","Maulisia_microlepis","Normichthys_operosus","Searsia_koefoedi","Sagamichthys_schnakenbecki" ) ~"Platytroctidae", species %in%c("Sigmops_bathyphilus","Gonostoma_elongatum") ~"Gonostomatidae", species %in%c("Argyropelecus_hemigymnus","Maurolicus_muelleri","Argyropelecus_olfersii" ) ~"Sternoptychidae", species =="Anoplogaster_cornuta"~"Anoplogastridae", species %in%c("Arctozenus_risso", "Paralepis_coregonoides") ~"Paralepididae", species =="Bathylagus_euryops"~"Bathylagidae", species =="Cyclothone"~"Gonostomatidae", species =="Derichthys_serpentinus"~"Derichthyidae", species =="Eurypharynx_pelecanoides"~"Eurypharyngidae", species =="Evermannella_balbo"~"Evermannellidae", species =="Lestidiops_sphyrenoides"~"Lestidiidae", species =="Melanostigma_atlanticum"~"Zoarcidae", species %in%c("Photostylus_pycnopterus","Xenodermichthys_copei") ~"Alepocephalidae", species =="Serrivomer_beanii"~"Serrivomeridae" ) )%>%mutate(order =case_when( family =="Myctophidae"~"Myctophiformes", family %in%c("Stomiidae","Gonostomatidae", "Sternoptychidae") ~"Stomiiformes", family %in%c("Platytroctidae","Alepocephalidae") ~"Alepocephaliformes", family =="Anoplogastridae"~"Trachichthyiformes", family %in%c("Paralepididae","Evermannellidae","Lestidiidae") ~"Aulopiformes", family =="Bathylagidae"~"Argentiniformes", family %in%c("Derichthyidae","Serrivomeridae") ~"Anguilliformes", family =="Eurypharyngidae"~"Saccopharyngiformes", family =="Zoarcidae"~"Perciformes", ) )
## Summary of the assemblages * species data.frame ----asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)asb_sp_fish_occ <- asb_sp_fish_summ$"asb_sp_occ"htmltools::tagList(DT::datatable(asb_sp_fish_occ))
2.2 Computing distances between species based on functional traits
We have non-continuous traits so we use the Gower distance(metric = “gower”) as this method allows traits weighting.
2.3 Building functional spaces and chosing the best one
2.3.1 Computing several multimensional functional spaces and assessing their quality
mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).
2.3.3 Testing the correlation between functional axes and traits
Code
sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"# As we have 26 traits we have to split the df to see correlation between functional axes and traits # first set ----fish_traits_1 <- fish_traits%>%select(1:9)fish_tr_faxes <- mFD::traits.faxes.cor(sp_tr = fish_traits_1, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
1 eye_size PC1 Linear Model r2 0.164 0.0087
2 eye_size PC2 Linear Model r2 0.174 0.0067
3 eye_size PC3 Linear Model r2 0.101 0.0425
4 eye_size PC4 Linear Model r2 0.162 0.0091
5 orbital_length PC1 Linear Model r2 0.406 0.0000
6 orbital_length PC2 Linear Model r2 0.104 0.0394
10 gill_outflow PC2 Linear Model r2 0.513 0.0000
13 oral_gape_surface PC1 Linear Model r2 0.126 0.0225
14 oral_gape_surface PC2 Linear Model r2 0.500 0.0000
18 oral_gape_shape PC2 Linear Model r2 0.150 0.0123
24 oral_gape_position PC4 Linear Model r2 0.228 0.0016
26 lower_jaw_length PC2 Linear Model r2 0.685 0.0000
29 head_length PC1 Linear Model r2 0.330 0.0001
30 head_length PC2 Linear Model r2 0.382 0.0000
34 body_depth PC2 Linear Model r2 0.353 0.0000
35 body_depth PC3 Linear Model r2 0.196 0.0038
Code
## Plot ----fish_tr_faxes$"tr_faxes_plot"
Code
# second set ----fish_traits_2 <- fish_traits%>%select(10:18)fish_tr_faxes_2 <- mFD::traits.faxes.cor(sp_tr = fish_traits_2, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
2 pectoral_fin_position PC2 Linear Model r2 0.278 0.0004
5 pectoral_fin_insertion PC1 Linear Model r2 0.282 0.0004
6 pectoral_fin_insertion PC2 Linear Model r2 0.406 0.0000
9 transversal_shape PC1 Linear Model r2 0.116 0.0297
11 transversal_shape PC3 Linear Model r2 0.218 0.0021
14 caudal_throttle_width PC2 Linear Model r2 0.401 0.0000
19 dorsal_fin_insertion PC3 Linear Model r2 0.190 0.0044
23 eye_position PC3 Linear Model r2 0.191 0.0043
26 operculum_volume PC2 Linear Model r2 0.104 0.0393
32 ventral_photophores PC4 Kruskal-Wallis eta2 0.590 0.0000
33 gland_head PC1 Kruskal-Wallis eta2 0.245 0.0011
Code
## Plot ----fish_tr_faxes_2$"tr_faxes_plot"
Code
# third set ----fish_traits_3 <- fish_traits%>%select(19:26)fish_tr_faxes_3 <- mFD::traits.faxes.cor(sp_tr = fish_traits_3, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value"<0.05), ]
Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units
null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurrence frequency)
FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.
FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.
FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.
FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.
calcul des indices avec package FD
correction distance matrix (obtenue avec Gower) = “none” or “lingoes” transformation avec “sqrt” ne permet pas de transformer la matrice en “euclidean”
randomisation des traits des espèces par couche de profondeur
interprétation : La richesse fonctionnelle (FRic) augmente avec la profondeur pour atteindre un maximum en bathypélagique, ce qui correspond à nos anciens résultats. Cependant, les valeurs de SES FRic sont également positives en épipélagique, ce qui indique une richesse fonctionnelle plus élevée qu’attendue au hasard. En revanche, les SES FDIS sont significativement plus faibles en épipélagique, ce qui indique que les espèces représentant le plus de biomasse sont similaires sur le plan fonctionnel.
comment faire pour tester la symétrie et la normalité des valeurs de SES, suelement sur les quatres valeurs que j’obtiens par indices ?
Code
# Model parameters ----# Number of simulationsn_simulations <-9# correction method to use when the distance matrix cannot be represented in a Euclidean spacecorr_method <-"lingoes"# "none" # Depth layersdepth_layers <-rownames(depth_fish_biomass)# List of indices to be calculatedindices <-c("FRic", "FDis", "FDiv", "FEve") # Matrices for storing observed and simulated resultsdbFD_result_obs <-matrix(NA, nrow =length(depth_layers), ncol =length(indices), dimnames =list(depth_layers, indices))dbFD_result_sim <-array(NA, dim =c(length(depth_layers), n_simulations, length(indices)), dimnames =list(depth_layers, paste0("Sim.", 1:n_simulations), indices))# Function to randomize traits while preserving factor levels----randomize_traits <-function(traits_mat) { randomized_mat <- traits_mat# Randomisation of columnsfor (trait incolnames(traits_mat)) { randomized_mat[, trait] <-sample(traits_mat[, trait], replace =FALSE) }return(randomized_mat)}# Loop on each depth layer ----for (layer in depth_layers) {# Select the species present in the selected depth layer species_in_layer <-colnames(depth_fish_biomass)[depth_fish_biomass[layer, ] >0]# Filter lines (traits) and remove constant columns traits_layer <- fish_traits[species_in_layer, , drop =FALSE] %>%select(where(~length(unique(.)) >1)) # If no variable lines remain, we move on to the next layerif (ncol(traits_layer) ==0) next# Correct formatting of biomass_layer with depth layer in rownames# Converts biomass data to matrix, which is required by FD::dbFD() biomass_layer <-matrix(as.numeric(depth_fish_biomass[layer, species_in_layer, drop =FALSE]), nrow =1, dimnames =list(layer, species_in_layer))# Calculation of observed indices with FD::dbFD dbFD_result <- FD::dbFD(x = traits_layer, w.abun =TRUE, m =4, a = biomass_layer, messages=F,corr = corr_method)# Storage of observed values dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv # Simulations for each layerfor (sim in1:n_simulations) { randomized_traits <-randomize_traits(traits_layer) dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun =TRUE, m =4, a = biomass_layer, messages=F,corr = corr_method)# Storage of simulated values dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv } }# Calculation of the SES (Standardized Effect Size) ----SES_results <-data.frame(depth_layer =rownames(dbFD_result_obs))for (index in indices) { meanNullFD <-rowMeans(dbFD_result_sim[, , index], na.rm =TRUE) sdNullFD <-apply(dbFD_result_sim[, , index], 1, sd, na.rm =TRUE) SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD}# Plot ----results_df <- tidyr::pivot_longer( SES_results,cols =-depth_layer,names_to ="index",values_to ="SES_FD")results_df$depth_layer <-factor( results_df$depth_layer,levels =c("Epipelagic","Upper mesopelagic","Lower mesopelagic","Bathypelagic" ))results_df$index <-factor( results_df$index,levels =c("FRic","FDis","FDiv","FEve"),labels =c("Functional richness","Functional dispersion","Functional divergence","Functional evenness" ))ggplot(results_df, aes(x = depth_layer, y = SES_FD, color = depth_layer)) +geom_point(size =3) +facet_wrap(~index) +geom_hline(yintercept =c(1.96, -1.96), linetype ="dashed", color ="gray40", size=0.8) +scale_color_manual(values =c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +labs(x ="", y ="Standard Effect Size (SES)") +theme_light() +theme(axis.text.x =element_blank(),strip.text.x =element_text(size =14, color ="black"),strip.background =element_rect(fill ="white"),axis.title =element_text(size =13),axis.text =element_text(size =13)) +guides(col="none", fill="none")
# test_norm$depth_layer <- factor(# test_norm$depth_layer,# levels = c(# "Epipelagic",# "Upper mesopelagic",# "Lower mesopelagic",# "Bathypelagic"# )# )# # # Density plot with normal density overlay# ggpubr::ggdensity(test_norm, x = "SES", fill = "gray", col = "white") +# facet_grid(depth_layer~index, scales = "free") +# ggpubr::stat_overlay_normal_density(color = "red", linetype = "dashed", size = 1) + # theme_light() +# labs(title = "Distribution of values of the random communites in compariaosn to the normal distribution",# x = "Values", # y = "Density")